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ABSTRACT 

We use the triaxial halo model formalism of lSmith fc Watt^ ^^200?^^ to investigate the 
^ , impact of dark matter halo shapes on the cosmological bispectrum. Analytic expres- 

sions for the dark matter distribution are derived and subsequently evaluated through 
numerical integration. Two models for the ellipsoidal halo profiles are considered: a 
. toy model designed to isolate the effects of halo shape on the clustering alone; and 

the more realistic model of Jing & Suto (2002). For equilateral A:-space triangles, we 
show that the predictions of the triaxial model are suppressed, relative to the spherical 
model, by up to ~ 7% and ^ 4% for the two profiles respectively. When one considers 
' the reduced bispectrum as a function of triangle configuration it is found to be highly 

, sensitive to halo shapes on small scales. The generic features of our predictions are 

■ that, relative to the spherical halo model, the signal is suppressed for fc- vector con- 

' figurations that are close to equilateral triangles and boosted for configurations that 

lO ■ are colinear. This appears to be a unique signature of halo triaxiality and potentially 

' provides a means for measuring halo shapes in forthcoming cosmic shear surveys. 

The galaxy bispectrum is also explored. Two models for the halo occupation dis- 
tribution (HOD) arc considered: the binomial distribution of Scoccimarro et al. (2001) 
I ' and the Poisson satellite model of Kravtsov et al. (2004). Our predictions show that 

O I the galaxy bispectrum is also sensitive to halo shapes, although relative to the mass 

> ' the effects are reduced. The HOD of Kravtsov et al. is found to be more sensitive. This 

^ I owes to the fact that the first moment of the occupation probability is a steeper func- 

• • , tion of mass in this model, and hence the high mass (more triaxial) haloes arc more 

^ ' strongly weighted. Interestingly, the functional form of the configuration dependent 

k>( I bispectrum is, modulo an amplitude shift, not strongly sensitive to the exact form 

. of the HOD, but is mainly determined by the halo shape. However, a combination of 

' measurements made on different scales and for different fc-space triangle configurations 

is sensitive to both halo shape and the HOD. 

Key v^rords: Cosmology: theory - large scale structure of Universe - Galaxies: grav- 
itational clustering 



1 INTRODUCTION 

In our current picture of structure formation, small Cold 
Dark Matter (CDM) density fluctuations, seeded during an 
early epoch of cosmic inflation, collapse through gravita- 
tional instability. From an initial spectrum of fluctuations 
that is almost scale-invariant, the collapse proceeds in a hi- 
erarchical way with fluctuations on small scales collapsing at 
early times to form dense clumps of dark matter (haloes), 
and with those on larger scales collapsing at later times. 
This gives rise to a web- like network of filam ents and clus- 
ters - the 'cosmic web ' jBond fc Mverslll99d) . Current ob- 
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servational constraints strongly suggest that the growth of 
structure takes place in a Universe that is flat and which 
at late times is dominated by an unknown Dark Energy 
component that dr ives an accelerated Universal expansion 
iSpergel et al .120031) . Building a complete statistical descrip- 
tion of the clustering pattern of the dark matter is conse- 
quently of great importance: this may provide detailed in- 
formation concerning the physics of the dark matter, the 
statistics of the primordial fluctuations and also valuable 
insight into the nature of dark energy. Direct comparison of 
the dark matter distribution with that of the galaxies reveals 
the salient physics of galaxy formation. 

However, the dynamics of structure formation are 
nonlinear and analytic solutions for the equations of 
motion are only tractable for highly idealized cases. 
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Robust calculation of the CDM density field therefore 
requires recourse to numerical integration via A''-body 
methods. Over the past decades significant effort has 
been invested in performing and analyzing large numerical 
simulations. The results of this have led to a number of 
semi-empirical formulae that describe the dark matter 
halo phenom enology to high accura, C Y. In particular, hal o 
abundances llSheth fc Tormcn ' l99gl: I Jen kins ct al.' '200ll) 
large -scale bias Hmo k White! Il99g~ S licth & Tormen 
12221), density str u ctures JNayarro. Frenk fc Whitf 
1997 - hereafter NFWt iBullock et alJ I2001I: iPower et al' 
20031: IJing fc Sutdl2002 - hereafter JS02ll . mte rnal velocity 
distributions and large scale velocity bias iSheth et alJ 
[2001a b ; Sheth & Diafcrio 200llh are now well characterized. 
More recently, attention has been focused on quanti- 
fying halo shapes (JS02), substructure phenomenology 
jMoore et al J 119991 : iKlvpin et al1ll999l: iGhigna et all 1200011 
and effects on the halo due to baryon physics Jcnedinet^Ll 
I2OO4, : Kaza ntzidis et al. 2004). Collectively, these ideas have 
lead to the development and evolution of the halo model 
paradigm for structure formation. 

In the most basic formulation of the halo model JSeliakI 
l200d:IPeacock fc S mith 2000; M a fc Frvll2000h . one assumes 
that all of the dark matter in the Universe is confined to 
virialized, spherical haloes that possess some distribution in 
mass, some universal density structure and a large-scale bias 
derivable from perturbation theory, all of which are charac- 
terizable through the halo mass. The clustering statistics 
over all scales then conveniently break down into terms in- 
volving correlations between haloes and correlations within 
haloes. This approach has had a high degree of success in 
modelling the low-order spatial clustering statistics of the 
mass. Most advantageously, the halo model can also be used 
to make predictions for the clustering statistics of the galaxy 
distribution, which in general can be thought of as a bi- 
ased sampling of the underlying mass distribution. Indeed, 
in the halo picture the concept of bias is replaced with 
that of the Halo Occupation Distribution (HOD), which 
dictates the nu mber and position of gala xies within dark 
matter haloes iBerlind fc Weinberg! 120031 . Moreover, the 
HOD connects in a satisfying manner to bot h the theoret- 
ical predictions of galaxy formation models (Benson et al. 
boOQ; S_cliak2000; ScQccim a rro et al. 200 1 -_hcrcMtcr__SSHJ; 
iBerlind fc Weinb erg |200 2|.~and others ) and direct observa- 
tional realitv llPea cock fc Smithll200d| : iBerUnd fc Weinberg 
l2002l:lYang. Mo fc van den Boschll200a . and others'). 

Whilst the agreement between the halo model and nu- 
merical simulations is in general good at the two-point level, 
significant d ifferences have been reported a t the three-point 
level rSSHJ: iFosalba. Pan fc Szapudillioosl) . Although a re- 
cent study bv ^WM^^^^^i^2004^ showed that the spherical 
halo model could be modified to provide a better fit to sim- 
ulation data, if a finite volume correction and a somewhat 
arbitrary spatial exclusion scale for the haloes were taken 
into account. However, discrepancies still remain discernible 
on both small and large scales. It is currently believed that 
the various inconsistencies between the halo model and N- 
body simulations are manifestations of the breakdown of key 
assumptions in the model. Thus, given the significant inter- 
est in this halo based approach, it is of crucial importance to 
understand the validity of all its approximations. In particu- 



lar, these issues must be addressed if we are to take seriously 
the halo model as a tool for precision cosmology. 

One possible erroneous assumption is that dark mat- 
ter haloes are spherical, or rather, that we may work with 
the spherical averages of the density profiles. However, it 



has long been known that the haloes 
simulations are more closely described 


found in numerical 
is triaxial ellipsoids 


jBarnes fc Efstathiou 1987: 


Frenk et al. 


1988: Warren et alJ 



Il992h . Recently, JS02 have shown, through high resolution 
A'^-body simulations, that the inclusion of halo shape infor- 
mation into the modelling of the density structure leads to 
improved fitting functions. With these ideas in mind, some 
key questions can now be asked: To what extent do the 
shapes of the dark matter haloes affect the clustering statis- 
tics of the mass distribution? Is there any observable affect 
on the galaxy clustering? 

In a recent paper, ISmith fc Watt j J2005l hereafter 
SW05) laid out the foundations for performing halo model 
calculations with triaxial haloes: the 'triaxial halo model'. 
They then went on to calculate the importance of halo 
shapes on the dark matter power spectrum, the Fourier 
transform of the two-point correlation function. They re- 
ported that the effect of halo triaxiality was to suppress the 
power spectrum on small/nonlinear scales at the level of be- 
tween ~ 5% and ~ 15%, dependent on the precise choice of 
density profile model. 

A number of other modifications to the halo model 
may lead to similar effects on the power spectrum. Partic- 
ular examples are the inclusion of a s tochastic concentra- 
tion parameter in halo density profiles JCoorav fc Hull200ll) 
and the treatment of halo su bstructure iSheth fc JairJl200; ; 
iDolnev. Jain fc Takadal2004 . More recently, some have sug- 
gested that changes to the halo boundary definition may 
affect significant changes to the clus tering iTakada fc JainI 
l2003al : lFosalba. Pan fc Szapudill2005^ . Clearly, if one wishes 
to disentangle these effects one must look to higher order 
statistics. In this paper we apply the triaxial halo model for- 
malism of SW05 to the task of computing the bispectrum, 
the Fourier space analog of the three-point correlation func- 
tion. The bispectrum is considered to be the lowest order 
spatial statistic that is sensitive to the shapes of structures. 
We will demonstrate that this measure can indeed be used 
to differentiate between changes to halo density profiles and 
changes to halo shapes. 

The paper breaks down as follows: In Section |5| we re- 
view some theoretical notions and lay down the bones of the 
triaxial halo model formalism. In Section |21 we provide the 
details of our derivation of the three-point correlation and 
bispectrum for the mass, presenting the results in Section 
|1| In Section 13 we switch from mass clustering to galaxy 
clustering and repeat the analysis for galaxies. Finally, we 
discuss our findings and draw our conclusions in Section |S] 

Throughout we assume a fiat, dark energy dominated 
cosmological model with equation of state manifest as a cos- 
mological constant. We take flm ~ 0.3 and ^Ide ~ 0.7, 
where Qm and Qde are the ratios of the density in dark mat- 
ter and dark energy to the critic al density, respectively. We 
use th e linear power spectrum of lEfstathiou. Bond fc White! 
jl99^ with normalization erg = 0.9 and shape parameter 
r = ^l,nh — 0.21, where h — 0.7 is the dimensionless Hubble 
parameter. 
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2 THEORETICAL BACKGROUND 

2.1 three-point spatial clustering statistics 

In this paper we are concerned with the three-point spa- 
tial statistics of the density contrast field S{r). This field is 
defined through the relation 



(1) 



where p(r) is the physical density and p is the density of the 
background. The three-point correlation function, (, is de- 
fined as the ensemble average of S measured at three points 
in space. This can be written 



C(ri,r2,r3) = (<5(ri)5(r2)5(r3)) 



(2) 



where the angled brackets denote the ensemble average. For 
homogeneous random fields ( obeys translational invariance: 



C(ri,r2,r3) = C(ri +ro,r2 +ro,r3 -fro) 



(3) 



where ro is an arbitrary position vector. For isotropic ran- 
dom fields, also obeys rotational invariance: 



C(ri,r2,r3) = C(TZri,TZr2,Ur3), 



(4) 



where TZ represents an arbitrary coordinate rotation. For 
universes that obey the cosmological principle these two 
properties must hold. Lastly, ( must also be invariant under 
parity transformations: 



C(ri,r2,r3) = C(-ri, -r2, -r3) 



(5) 



In what follows, we will not deal explicitly with (, but 
instead work with its Fourier transformed counterpart, the 
bispectrum B, 

{5(ki)5(k2)<S(k3)> = (27r)'B(ki,k2,k3)<S''(ki+k2+k3),(6) 

where (5(k) is the Fourier transform of 5(r) and is the 
Dirac delta function. From the argument of the delta func- 
tion we obtain the triangle condition, i.e that the sum of the 
three wave- vectors ki, k2 and k3 form the null vector. The 
bispectrum and three-point correlation function are explic- 
itly related through 

dki dk2 



C(ri,r2,r3) 



B(ki,k2,-ki -k2 



(27r)3 (27r)3 

X exp [-iki ■ (ri - r3)] exp [-ik2 ■ (r2 - rs)]. (7) 
From the reality of i^, we find that 

B(ki,k2,k3) =B*(-ki,-k2,-k3) , (8) 

where * denotes complex conjugation and henceforth the tri- 
angle condition is to be assumed. If the 5-field is statistically 
homogeneous and isotropic, then following equation the 
bispectrum also obeys parity invariance: 



B(ki, k2, k3) = B(-ki, -k2, -k3 



(9) 



This, combined with the reality constraint of equation © 
implies that the bispectrum is real. Lastly, following equa- 
tion 1^ the bispectrum is also rotationally invariant: 



B(ki,k2,k3) = B(7^kl,7^k2,7^k3 



2.2 Ellipsoidal systems 



(10) 



In what follows we will be concerned with the cluster- 
ing properties of triaxial dark matter haloes. We there- 



fore define some imp ortant relations for these objects (see 
IChandrasekharl Il969l for a full treatise). Consider a het- 
erogeneous triaxial ellipsoid that has semi-axis lengths a, 
b and c, where a < b < c, and orthogonal principle axis 
vectors (ea,ei,,ec). We define the triaxial coordinate sys- 
tem (7?, 0, $) with respect to the principle axes of the halo, 
where Be is taken to be in the z-direction of a standard 
Cartesian system. The radial parameter R traces out thin 
iso-density shells, or homoeoids, and the parameters B and 
$ are the polar and azimuthal angles respectively. In this 
system the Cartesian components can be written 

X = — J? cos $ sin O ; y = -i?sin $ sin Q ; z — J? cos O. (11) 

It is to be noted that the ellipsoidal angles differ from those 
of the spherical coordinate system. The parameter R can be 
related to the Cartesian coordinates and axis ratios through 

2 ,2 



R 



y 



(12) 



~^ b'^ ~^ c? 

The benefits of this choice of coordinate system are now 
apparent: if all homoeoidal shells are concentric, and if one 
picks coordinates with the same axis ratios as the triaxial el- 
lipsoid, then the density run of the ellipsoid can be described 
by a single parameter: 



p(r) ^ piji) 



(13) 



Furthermore, the mass enclosed within some iso-density cut- 
off scale -Rcut, can be obtained most simply by 



ab 



dRR^p{R) , 



(14) 



M = I drp(r) 

J flout (r) 

where the ellipsoidal coordinates have allowed us to circum- 
vent the problem of evaluating complicated halo boundaries. 

It is now also convenient to define what we mean by 
a halo: any object that has a volume averaged over-density 
200 times the background density is considered to be a grav- 
itationally bound halo of dark matter. This leads directly to 
the following relation between the mass, radius and axis ra- 
tios 



M2, 



4 p3 /afc 



(15) 



The above definition was adopted in order t o be consistent 
with the mass- function of ISheth fc TormenI (^9^. 



2.3 Triaxial halo model 

We next summarize the theoretical ingredients of the triax- 
ial halo model presented in SW05. Consider a density field 
comprised entirely of ellipsoidal dark matter haloes of dif- 
ferent masses. The density run of each halo is specified by 
some universal profile, the exact details of which are not 
important at this stage save that the general properties dis- 
cussed in the previous section are satisfied. Each halo may 
therefore be characterized by a set of stochastic variables. 
These are the position vector for the halo centre of mass, 
X, the mass M, the principle axis lengths (a,fe, c), and the 
direction vectors for these axes (e^, et,, ec). We express these 
last two sets of variables compactly by using the notation 
a = (a,6, c) and £ = {ea,eb,ec). Thus, the density at any 
point r can be expressed as a sum over the A'^ haloes that 
comprise the field 
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(16) 



where U is the mass normalized density profile and where 
sub-scripts i denote the characteristics of the i"' halo. 

The correlation functions of the density field follow di- 
rectly from taking ensemble averages of products of the den- 
sity at different points in space. To compute the ensemble 
averages we integrate over the joint probability density func- 
tion (hereafter PDF) for the haloes that form the density 
field, and sum over the probabiliti es for obtaining the A'^ 
haloes fsee iMcClelland fc Silklll977l for a similar approach) 



(...) 



p{N,\V) 



J Y[ dMi dx.i dsLi d£i 



X p(Mi, .., AfjVj. ,xi, ..,xjVj ,ai, .., slnj , £i, £ \Nj) .(17) 

Provided the volume of space considered is large, then 
p{N\V) is very sharply spiked around N = nV 1, where n 
is the mean number density of haloes. We restrict our study 
to this case only. 

The integrals over £ in equation 11711 represent averages 
over all possible orientations of the halo. The orientation of 
the halo frame can be specified relative to a fixed Cartesian 
basis set through the Euler angles. These represent succes- 
sive rotations of the halo frame about the ^:-axis by a, the 
jy'-ax is by /3, and the 2"-axis by 7 (see lMathews fc Walked 
ll97(J) . Hence, the process of averaging over all halo orienta- 
tions can be performed by integrating over all Euler angles. 
We thus make the following transformation £ —> {a,l3,"/). 

Finally, we will require explicit expressions for the joint 
PDF for a single halo's characteristics, p(x, Af, f ,a). If we 
assume that a halo's orientation, position and mass are in- 
dependent random variables, and that the halo axis ratios 
are dependent on mass only (see e.g JS02), then we may 
write 

p(x,M,f,a) = i^p(£)p(a|M), (18) 
V n 

where n{M) is the mass function of haloes. If the orientation 
of a halo is uniformly random on the sphere then p{£) = 
I/Stt^ (see SW05 for further details). 



3 CLUSTERING OF TRIAXIAL HALOES 

3.1 The three-point correlation function 

We are now in a position to compute the three-point correla- 
tion function ("(ri,r2,r3) in the triaxial halo model. This is 
done by applying the machinery of equations ll(ill and 11711 . 
to the definitions given by equation (|5J in conjunction with 
equation Q - see also SW05. We find that, as for the case 
of the standard halo model, ^ can be expressed as the sum 
of three terms 



(19) 



where we have suppressed the position vector dependence 
for simplicity. Here C^^^, ^'^^ and C^^^ give the contributions 
for the case where: all three points are within the same halo; 
two points are within the same halo and the third is within 
a separate halo; all three points are within distinct haloes. 



Throughout, we refer to these contributions as the 1-, 2- and 
3-Halo terms. These have the following forms: 



^IH _ 1 



p387r2 



dM dx da d£ n(M)p(a|A/) 



]^{[/(M,x-r„a,£)} 



(20) 



^2H _ 



1 



JJ^ {dMi dx.i dsLi d£i n{Mi) p(a.i\Mi 



p3(87r2)2 

X U(xi~ri,Mi,ai,£i)}Mi M2 
X [/(xi -r3,A'/i,ai,£:i) f (1,2) -f eye 

1 



(21) 



]^{dAf, dxi dai d£i Min{Mi)p{ai\Mi) 



X C/(x, -r„A'/„a„£:0} C'(l,2,3) 



(22) 



In deriving this result we have used the fact that the joint 
PDF for the characteristics of three haloes may be rewritten 



p(l, 2, 3) = p(l) p(2) p(3) [1 + r (1, 2) + r (2, 3) 
+r(3,l) + r(l,2,3)] , 



(23) 



where we have used the short-hand notation p(l) = 
p(xi, Afi, ai, f 1) etc. The quantity p(l) is the 1-Halo PDF, 
as given by equation 11811 . while ^^(1,2) and ("^(1,2,3) are 
the two- and three-point seed correlation functions for tri- 
axial haloes. 

In the work of S W05 it was shown that the effects of halo 
alignments on the mass clustering statistics are insignificant, 
even for the highly symmetrized case where all of the haloes 
are perfectly aligned. We therefore take the orientation and 
shape of each halo to be statistically independent of each 
and every other halo. Under this condition the joint halo 
PDF for the three-point halo properties reduces to 

p(l,2,3) =p(l)p(2)p(3)[l + r(Mi,A^2,xi,X2) 

( Af2 , Af3 , X2 , X3 ) 4- f ( , A'/iX3 , Xi ) 
-hC''(Ml,A^2,A'/3,Xl,X2,X3)] , (24) 

where the seed correlation functions can now be specified in 
the usual way iIMo fc Whitelll99fil: Isheth fc Tormenlll99fll^ . 



3.2 The Bispectrum 

The bispectrum is most directly obtained through Fourier 
transforming equation l|19|l . From the linearity of C,, we see 
that B is also composed of three terms 



(25) 



and where we have adopted the short-hand notation -B123 = 
_B(ki, k2, k3). In the subsequent sub-sections we derive ex- 
plicit relations for these quantities. 

3.2.1 The 1-Halo term 

Fourier transforming equation 12UII leads us to the 1-Halo 
contribution to the bispectrum. This is compactly written 
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= 4 / dMA'I^n{M)Wi23{M) , 



(26) 



where we have introduced the (three-point) window function 
Wi23{M) = ^ I dad£p{a\M) J| (7(k„ M, a, f), (27) 

•' i={l,2,3} 

where U(ki, M, a,£) is the Fourier transform of the mass 
normalized triaxial halo profile. In its present form equa- 
tion 12711 is of limited practical use, since its computation 
requires one to solve a very high dimensional numerical in- 
tegral. In Section 13.31 we show how this expression can be 
greatly simplified. 

3.2.2 The 2-Halo term 

The Fourier transform of equation 12H gives the 2-Halo con- 
tribution to the bispectrum. This is succinctly written 

b!23 = ^^^J n {dM,n{M,)biAh)} 



»={1,2} 

X Mf Ma Wi3(Mi) WiiMi) + eye 



(28) 



where Wi and Wij axe the 1- and two-point window func- 
tions defined by the relations 

W,{M) = J dap(a|M) J d£U{ki, M, &,£) ; (29) 
and 

Wij{M) = ^ j dad£p{a\A4) J| f/(k„ M, a, f) ; (30) 

and where we have written the Fourier transform of the 
halo seed two-point correlation function as P"{k, Mi, M2) = 
6i(Afi) fei(M2) -F^(fe), where -f^ is the linear power spec- 
trum and 61 (M) is the first linear halo bias parameter 
JMo fc WhitellQQdISheth fc Tormenll99g^ . Again, these ex- 
pressions will be simplified following the discussion in Sec- 
tion |^] 



3.2.3 The 3-Halo term 

Finally, the 3-Halo contribution to the bispectrum is ob- 
tained through Fourier transforming equation 122L 

1 r ''^ 

Big = — \{{dM,M,n{M,)W,{Mi)} 



X Bl23(Mi,M2,Mi) 



(31) 



where Wi{M) is given by equation I29II . The function 
Bl2i{M\, M2, M3) is the bispectrum of halo seeds and is 
given by (SSHJ) 



B!23(Ml,M2,M3) = 61 (Ml) 61 (M2) &l(M3) 



i>123 + 



61 (M3) 



Phiki) Phiki) + eye 



, (32) 



where &2(M) is the second order halo bias factor. The quan- 
tity -B123 is the second order Eulerian perturbat ion theory 
bispectrum llFrvlll984l : Ijain fc BertschingeJll994) . 



Bm = 2F2(ki,k2)PL(fci)PL(fc2) + eye , (33) 
and where 

F2(ki,k2) = ^ + icos6'i2(fci/fc2 + fc2/fci) + ^eos^0i2 ,(34) 
with cos 612 — ki ■ k2/fcifc2 . 



3.3 Simplification of halo window functions 

Having written down the basic results for the bispectrum, 
we now consider in more detail the expressions for the halo 
window functions defined by equations 12711 . 1291 and 13UII . 
So far, these have been left in terms of U(ki, M,a,£), the 
Fourier transform of the triaxial profile. As mentioned, nu- 
merical evaluation of these quantities in their present form 
would be very cumbersome, owing to the dependence of the 
halo profile on the orientation and shape of the halo. This 
means that the (full 3-D) Fourier transforms cannot be pre- 
computed but instead must be evaluated inside a 5-D in- 
tegral. The following arguments serve to simplify these ex- 
pressions considerably. 

Consider first the one-point window function given by 
equation 12911 . on inverse Fourier transforming the halo pro- 
file we find. 



X exp(— iki ■ r) . 



d£ dr U{r, M, a, £) 



(35) 



If we now construct a Cartesian coordinate system about the 
direction vectors of the semi-axes of the halo, then we may 
transform from that system to a system of ellipsoidal coor- 
dinates where the isodensity surfaces of the halo are given 
by equation 1121 . Thus, in this new basis set the density 
profile simply becomes a function of the ellipsoidal radial 
parameter R: U{r) —> U{R). 

The orientation average can now be performed di- 
rectly by rotating the reference Cartesian coordinate sys- 
tem through all possible values of the Euler angles. Hence, 
the components of the fc-vector are modified through each 
infinitesimal rotation {5a, 5 cos (3, 5"/), which leads us to 



m(M) 



dap (a I A/) 



ab 



dRR^U{R,M, a) 



X J dad{cosf3)dj J dRexp [-ik| ■ r(R)] , (36) 

where k' = TZ{a, /3, 7)k and TZ{a, (5, 7) is the rotation matrix 
(see Appendix^for an explicit definition). Following SW05 
the integral over R can be performed analytically using the 
relation 



dRexp [-ik • r(R)] = 47rjo [kRfiOk, < 



(37) 



where jo is the zeroth order spherical Bessel function and 
f^{dk,4'k) = cos^ 9k +sm^ 9k ( ^ cos^ (j)k + ^ sin^ <j)k] ■ (38) 



Finally, on substituting equation 1371 into 1361 we find 
m(M) = ^ [ dap(a|A/)^ f dad(cos/?)d7 , 
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Figure 1. The triaxial halo window functions compared with the 
equivalent spherical halo window functions. In the top panel, thick 
lines represent triaxial haloes and thin lines spherical. Different 
line styles represent Wi, W12 and W123 as defined by equations 
I39i . I4.'il and 1441 . The lower panel presents the ratio of the 
triaxial window functions to the spherical. The line styles are as 
in top panel. 



X J dRR^U{R,M)jo[hRfiei^A'k,)] 
where we have from the rotation matrix that 



— arccos 



Dfc = arccos 



k 



V 



k sin 6', 



(39) 



(40) 



(41) 



T] = {C fiC ctC'^ SctS^^kx'\~{C fiSaC^'\'CaS^^ky iS'/^C^-^/b^ (42) 

and where we have adopted the short-hand notation Cx = 
cos a:: and Sx = sinx. Importantly, ifa/c = 6/c=l, then the 
window function simplifies to that of the standard spherical 
halo model. 

Similarly, applying this procedure to the functions Wij 
and VF123, we find: 



Wij{M) 



(470 
87r2 



dap(a|M) [ da d{cos jS) d-y 



X Yl I dR.RlU{R., M, a) }o[k.R.f{9'^,^,(t>'k jy, (43) 
and 



(4^ 
87r2 



dap(a|M) (^) 



da d(cos P) d7 



n 



dR,Rl U{R,,M, a)io[k,R.f{e'^^ , J].(44) 



In accordance with the set of conditions on the bispectrum, 
presented in Section [2. II we see that the window functions, 
specified by equations 139L 14311 and 14411 . are invariant un- 
der arbitrary rotations of the k vector triple, and that they 
are also real functions of k. 

As a point of practice for computing these functions, 
it is necessary to specify the initial coordinates of the k- 
vector triple. There is no loss of generality here, owing to 
the orientation average of the halo. We therefore choose ki 
to be aligned with the z-axis of the reference coordinate 
system. The orientation of k2 is then specified in spherical 
polar coordinates so that the polar angle gives the angle 
between ki and k2. The azimuthal angle of k2 is then, from 
equations 14UII and 1421 . entirely arbitrary. Once k2 has been 
specified, ka is then fixed through the triangle constraint. 
The window functions, and subsequently the bispectrum, 
then simply become functions of |ki|, |k2| and the angle 
9 = arccos (ki • k2/fcifc2), as required for homogeneous and 
isotropic random fields. 

Figure shows a calculation of the window functions 
given by equations I39L I43II and I44II for equilateral trian- 
gle configurations, i.e |ki| = |k2| = k and 612 — 27r/3. In 
this calculation we have taken the ellipsoidal density pro- 
file to be the triaxial NFW model (see Section 14.111 . and 
for illustration we have taken the joint axis ratio PDF to 
be simply the product of two delta functions: p{a/c, b/c) = 
5^{a/c — 0.5) 3^{b/c — 0.5). This constrains all haloes to take 
the form of prolate ellipsoids. The top panel shows the two 
usual characteristics: on large scales the window function ap- 
proaches unity and on small scales it rapidly decays with os- 
cillatory features present. To show the effect of halo triaxial- 
ity on these window functions in more detail, we compute the 
ratio of these quantities with the windows for the spherical 
model that result when p{a/c, b/c) = 5^(a/c— 1) 5^{b/c~l). 
The results of this are presented in the lower panel of the 
same figure. Clearly, the function VK123 is much more sen- 
sitive to the shape of the haloes than either Wi or W12. 
This leads us to conclude that the 1-Halo term becomes in- 
creasingly sensitive to the halo shape as the order of the 
clustering statistic increases. 



4 RESULTS: MASS BISPECTRUM 

In this section we evaluate the bispectrum given by the sum- 
mation of equations 12611 . 12811 and 1311 1. using our simpli- 
fied expressions for the window functions given by equations 
^ and (El. 



4.1 Computational details 

We investigate two models for the halo profile U{R, M,a.): 
the first is a toy model that allows us to explore how mod- 
ification of the halo shape alone affects the clustering - we 
refer to this as the 'triaxial NFW' model; the second is the 
more realistic density profile model of JS02, in which halo 
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Figure 2. The equilateral bispectrum. Top panel presents the dimensionless bispectrum as defined in equation 1461 . Thick lines represent 
the triaxial haloes and thin lines the equivalent spherical haloes. The dot, dot-dash and dash lines represent the 1-, 2- and 3-Halo terms 
respectively and the solid lines represent the total. Bottom panel presents the ratio of the triaxial halo bispectrum to the spherical. The 
line styles are preserved, but note that the dotted line represents the ratio of the triaxial 1-Halo term to the spherical halo model 1-Halo 
term, similarly for the 2- and 3-Halo terms. 



shape and central densities are not independent - we refer to 
this as the 'JS02 model'. For a more thorough explanation 
of these we refer the reader to SW05 and JS02, respectively; 
note that in SW05 the triaxial NFW model is referred to as 
the 'continuity model'. Common to both of these models is 
the density run, which is assumed to follow the NFW form 

^2 ' y^R/RoiM), (45) 

where R is the ellipsoidal radial parameter from equation 
112II . (5*'^""'(M) is the characteristic density of the halo and 
Ro{M) is the scale radius; in the JS02 model these are 
treated as independent variables. 

The PDF for the axis ratios, p( a|M), is taken from JS02 . 
We point out that in a recent paper. lLee. Jing fc Sutol (|200^ 
valiantly attempted to derive the halo axis ratio PDF from 
the peak theory of Gaussian random fields. However, their 
predictions appear to be in conflict with the results from 
A'^-body simulations and so we do not consi der their model 
further. We use the lSheth fc TormenI mass function 

and halo bias functions, and correct for the different defini- 
tions of halo mass as discussed in SW05. In order to focus 
purely on the effects of triaxial modelling, we do not draw 
the concentration parameter from a probability distribution, 
as advocated by JS02, but treat it as a deterministic variable 
related to the halo mass. 



Although in the previous section we greatly simplified 
our expressions for the halo window functions, we are still 
faced with the challenge of performing 7-D numerical inte- 
grals to compute the bispectrum. Evaluation of these ex- 
pressions through a set of serial quadratures is unfeasible on 
a standard workstation. Instead, we advocate the use of an 
efficient multi-dimensional algorithm, s uch as the Korobov- 
Conroy (fKorobov"l963'; 'Conrov* [ l967i) or Sag-Szekeres al- 
gorithm ( Sag & Szckcrcs 1964) . However, the integrals over 
the ellipsoidal radial parameter R require more care, due to 
the oscillatory nature of the integrand. To evaluate these 
we therefore employ a 1-D adaptive routine. Thus, for ex- 
ample, to evaluate the 1-Halo contribution to the bispec- 
trum for a given fc-space triangle (equation 1261 , we com- 
pute a 6-D numerical integral of a function that is itself the 
product of three 1-D integrals. Such numerical computations 
are achievable on a modern PC in a reasonable time frame. 
In Section [4.61 we present an independent determination of 
the triaxial halo model bispectrum, obtained through direct 
measurement from synthetic halo simulations. This provides 
a robust test of the veracity of the results produced by the 
multi-dimensional integrators. 
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Figure 3. Reduced bispectrum Q as a function of cosS = ki ■ k2/A;ifc2, where Q is defined through equation I47i and where we have 
considered the particular triangle configuration \k2\ = 2|fci|. The six panels show how the Q-9 curves evolve from large linear scales to 
small nonlinear scales; the sequence is shown for ^2 = (0.2, 1.0, 4.0, 10.0) hMpc~^ . Thick lines show the predictions of the triaxial halo 
model, thin lines are those for spherical haloes. 



4.2 Equilateral triangles 

Figure |5] (top panels) shows the predictions for the dimen- 
sionless bispectrum for equilateral triangles and for the tri- 
axial NFW and JS02 models. This is defined as 



AEQ(fe) 



47r 



■fcV-BEQ(fc) 



(46) 



where BEQ{k) is the bispectrum for equilateral triangles. 
In addition to results for triaxial haloes (thick solid lines), 
we also show the results for spherical haloes (thin lines). 
The spherical halo model predictions are obtained by set- 
ting the axis ratio PDF to be a product of delta func- 
tions as discussed in Section 13.31 It is evident that devia- 
tions from the spherical halo model appear to be very small 
for both profiles. In the bottom panels of Figure |2| we il- 
lustrate the differences more clearly by plotting the ratio 
A^q (triaxial) /Acq (spherical) for the 1-, 2- and 3-Halo term, 
and their sum, respectively. As may be expected from the 



window functions, it is clear that the biggest departure from 
the spherical model occurs in the 1-Halo term, which is at 
most ^ 7%. The 2- and 3-Halo terms also deviate by several 
percent, though only on scales below which their contribu- 
tion to the total bispectrum is insignificant. On larger scales 
there is no appreciable difference from the spherical model. 
Thus, for the purposes of halo model calculations that wish 
to measure the equilateral triangle bispectrum, it is a good 
approximation to neglect the halo shape information in the 
2- and 3-Halo terms. However, if one wishes to be accurate 
to the level of ~ 5% one must include shape information on 
the 1-Halo term. 

On the other hand, it is now clear that the discrepan- 
cies (~ 50%) between the measurements of the equilateral 
triangle bispectrum from numerical simulations and the pre- 
dictions for the halo model, no ted in the work of SSHJ and 
iFosalba. Pan fc Szapudil (j200^, can not solely be attributed 
to the break down of the spherical halo approximation. 



© 0000 RAS, MNRAS 000, 000-000 



The impact of halo shapes on the bispectrum in cosmology 9 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 

d/n d/n 



Figure 4. Dependence of Q{0) on the ratio k2/ki. Thick lines represent the triaxial halo model predictions and thin lines denote the 
spherical model. Solid, dash and dot-dash lines represent the cases where the triangle configurations: k2 = ki; k2 = 2ki; k2 = 3ki. 



4.3 Variation with triangle configuration 

We next consider the dependence of the bispectrum on the 
configuration of Fourier space triangles. This is achieved by 
fixing the lengths of ki and k2 and then varying the angle 
between them. Since the effects caused by halo triaxiality 
are relatively small, we choose to examine the more sensitive 
reduced bispectrum, or hierarchical amplitude: 

n,„ = ^ (47) 

^ ' ~ P{k^)P(k2) + P{k2)P{kz) + P{k3)P{kiy ^ ' 

where, for the case of the triaxial model predictions, P(k) is 
the triaxial halo model power spectrum (SW05). In Figure 
|H]we show the results for the triaxial NFW profile model on 
various scales. In all four panels we have taken one of the 
fixed sides of the triangle to be twice the length of the other 
fixed side: k2 = 2ki. The four panels show the results for 
the case where k2 = (0.2, 1.0,4.0, 10.0) /iMpc"^ 

The generic features of these curves have been examined 
in great detail in previous studies (Scoccimarro et al. 1998; 
SSHJ; Hou et al. 2005; Fosalba, Pan & Szapudi 2005). On 
large scales, it has been suggested that the variation in the 
configuration dependence can be attributed to coherent fil- 
amentary structures in the density field enhancing the con- 
tribu tion to Q from coUinear triangles (..Scoccim arro et alJ 
^9^. On intermediate scales, the density field begins to be 
dominated by collapsing objects that have almost isotropic 
spatial distributions, thus Q shows much reduced variation 
across the configuration. Finally, on small scales, Q is com- 
pletely determined by the structure and shapes of the dark 
matter haloes. This overall evolution is well reproduced in 
our Figure H 

Considering the relative differences between the triax- 
ial and spherical predictions in Fig. |3 we make the follow- 
ing observations: on large scales (fc2 = 0.2/iMpc~^; top left 
panel), where Q is dominated by the 3-Halo term, the af- 
fect of halo triaxiality is insignificant. This can easily be 



understood through considering the k —> limit of equa- 
tion 1391 1: On intermediate scales (fe — 1 /iMpc^^; top right 
panel), the magnitudes of the 1-, 2- and 3-Halo terms are 
comparable and one sees that there is a very small relative 
difference between the spherical and triaxial halo model pre- 
dictions. This difference originates from the growing 1-Halo 
term, there being no perceivable variations in the 2- or 3- 
Halo terms; On smaller scales (lower panels), Q is completely 
dominated by the 1-Halo term, and we see that there is a 
significant difference between the spherical and triaxial halo 
model predictions. The variation in Q across the set of tri- 
angle configurations, for the triaxial haloes, appears to be 
much flatter than the variation apparent in the spherical 
halo model. It should also be noted that the positions of the 
maxima and minima of Q are shifted. 

In Figure 2] we investigate these apparent differences 
between the spherical and triaxial halo model in closer de- 
tail. Fixing fei = 10/iMpc~^, we consider the cases where 
fe/fci = (1,2,3). Further, since the 1-Halo term dominates 
these scales, we focus in on the region over which it varies 
and ignore all other terms. In all cases we see that the tri- 
axial model predictions for Q, for both density profile mod- 
els, are suppressed with respect to the spherical model on 
angular scales where the fc-vectors are close to equilateral 
triangles, and amplified on scales where the fc-vectors are 
colinear and nonvanishing. Most importantly, from the fig- 
ures it is clear that the overall functional forms exhibited by 
the predictions depend little on the density profile model, 
modulo an amplitude shift, and that it is the shape of the 
haloes that dictates the shape of the curves. 

In the work of SSHJ, it was demonstrated that on 
small scales the predictions for the spherical halo model 
for Q{9) exhibited a convex configuration dependence, 
whereas the measurements from N-body simulations showed 
a much flatter dependence with 9. Indeed, SSHJ spec- 
ulated that this disagreement was likely the result of a 
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Figure 6. Dependence of Q{6) on halo mass. Again, thick Hnes represent triaxial model predictions thin lines spherical. 



breakdown in the assumption of halo sphericity. More re- 
cently these discrepancies have been highlighted by the 
work of iFosalba. Pan fc Szaptidil ll200?j ) , who demonstrated 
that even with a somewhat ad-hoc modification of the halo 
boundaries, the configuration dependence on small scales 
can not be replicated in the spherical halo model. On con- 
trasting these results with those above, we are lead to believe 
that the apparent discrepancies between theory and numer- 
ical simulation are likely attributed, in the main, to a break 
down in the spherical halo approximation. Moreover, it is 
highly unlikely that the exhibited configuration behaviour 
could be reproduced by other variations in the halo mod- 
elling prescription, e.g. through the stochasticity of the halo 
concentration parameter or through the inclusion of sub- 
structures. Such changes may alter the amplitude of Q{9), 
as between the left and right panels of Figure 3] but they 
should not change the overall shape. 



4.4 Dependence on halo shape 

We now explore in closer detail how sensitive the triax- 
ial halo model predictions are to the precise choice of the 
PDF for the axis ratios. This we do most directly by taking 
p{sl\M) as a product of Dirac delta functions (see Section 
13.31 . The sequence of curves in Figure |K| shows, in order of 
increasing thickness, the predictions for prolate dark matter 
haloes with axis ratios in the range 0.4 < a/c = b/c < 1.0, 
incremented in steps of a/c = b/c — 0.05. Again we have set 
fe/fci = 2 and fci = lQ/iMpc~\ For both the triaxial NFW 
profile and the JS02 model a generic trend is apparent; for 
spherical haloes, the maximum in Q{9) occurs for 9 = 2tt/3, 
with minima at 6^ = and n. As haloes become more el- 
lipsoidal the maxima in Q{9) is strongly suppressed while 
the minima are enhanced. For the most prolate objects the 
maxima in Q{9) lie at 6* = and vr, while a deep minimum 
forms between 9 = n/3 and 27r/3. 
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These trends can readily be understood by considering 
all possible placements of the fc-space triangle within the 
Fourier transformed halo profile ?7(ki, M, a, £). Varying the 
configuration of the triangle, one sees that the largest con- 
tribution to the bispectrum will come from triangles that 
most optimally fill the halo volume in fc-space. For a triax- 
ial ellipsoid, its Fourier transform is also ellipsoidally sym- 
metric. Thus, the 'optimal' triangle configuration is the one 
that matches most closely the halo's symmetry. For spherical 
haloes this clearly means equilateral triangles (or close to); 
whereas for very ellipsoidal objects (e.g. the thickest lines in 
Fig-UJi the optimum configurations shift to those triangles 
where the fc-vectors are coaligned. This explanation applies 
equally well to both the NFW and JS02 profiles since, by 
construction, they both share the same axis ratio PDFs. The 
apparent amplitude shifts in the JS02 predictions for Q{0) 
are due to the fact that in this model, both the central den- 
sity and concentration parameter depend on the halo axis 
ratios. 

4.5 Dependence on halo mass 

In the numerical simulation work of JS02 it was found that, 
for CDM models, the axis ratio PDF was conditional upon 
halo mass, and that more massive haloes were on average 
more triaxial than lower mass haloes. This can be under- 
stood from halo formation histories: high mass haloes form 
at late times in the Universe, whereas low mass haloes form 
at early times. Therefore, lower mass haloes have a longer 
period of time for particle orbits to undergo relaxation and 
become isotropic. Furthermore, since the haloes are assem- 
bled through the anisotropic accretion of lower mass haloes 
flowing along filaments, the halo shape at the time of forma- 
tion will be more closely described as a triaxial ellipsoid. It 
is therefore of interest to consider how the bispectrum pre- 
dictions depend on the abundances and masses of haloes. 

Figure |S] shows the contribution to Q{0) arising from 
four coarse bins in halo mass, for both the triaxial NFW and 
JS02 density profile models. Once again we set ^2/^1 — 2 
and = 20/1 Mpc~^. We take P{k) in equation 14711 as 
the power due to all halo masses, so that the sum of the 
contributions from the four bins equals the total Q. The 
main contribution to Q comes from haloes with masses in the 
range lO^*/i"^M0 < M < We see that both 

the triaxial and spherical model predictions exhibit a 'wave- 
like' form with 9. These waves are of similar amplitudes, 
but exhibit different 'phases'. The maximum amplitude for 
the triaxial halo model predictions occurs at = tt, with a 
second peak of slightly lower amplitude at 9 = 0, whereas 
for the spherical model, the peak amplitude is located at 
9 — 2n/3. These differences can be understood from the 
discussion surrounding Figure |S] 

The next most significant contribution comes from 
haloes with masses in the range 

lO"M0/i-\ These curves show less of the wave-like form 
exhibited in the higher mass haloes. Crucially, however, for 
the case of the spherical haloes the maxima and minima in 
Q occur at the same position as for the higher mass bin. For 
the triaxial predictions this is not the case, and the form 
of the curve has become almost consistent with that of the 
spherical halo model. This can be understood through the 
shifting of the mean of the JS02 halo axis ratio PDF to more 
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Figure 7. Test of the robustness of predictions of the triangle con- 
figuration dependent bispectrum. The star points show the mean 
bispectrum measured from an ensemble of 20 synthetic triaxial 
halo density fields ratioed to the predictions from the analytic 
spherical halo model. The circular points show the same, but for 
synthetic spherical halo density fields. The solid line shows the 
analytic predictions for the triaxial halo model. 

spherical haloes as the halo mass decreases. Finally, one can 
now see how for the triaxial halo predictions the summation 
of the contributions to Q from high and low mass haloes add 
de-constructively to give a relatively flat function across 9 
(cf. Fig.OJ. 

4.6 Robustness of predictions 

Owing to the relative obscurity of the high-dimensional in- 
tegrators that we have used to solve the integrals in Section 
13.21 we feel that it is necessary to demonstrate the robust- 
ness of their predictions. This we accompli sh through the use 
of syn theti c halo simulations. Followin g iPeacock fc Smithl 
J2000l) and lScoccimarro fc Shetl] i2002h . we generate direct 
realizations of the triaxial halo density held through popu- 
lating cubical regions of space with haloes that are sampled 
from the mass function, down to some limiting mass. Within 
each halo we then distribute equal mass particles in accor- 
dance with the required density structure, until the halo 
mass is reached. The spatial positions of the halo centres 
are chosen at random so that there is no large scale cluster- 
ing of halos. From this particle distribution we then directly 
measure the bispectrum. 

Figure |7| shows the measurement of the triaxial halo 
model and spherical halo model bispectrum in an ensem- 
ble of twenty synthetic simulations, ratioed with the ana- 
lytic predictions of the spherical halo model. These are then 
compared with the numerically integrated predictions of our 
analytic results for the bispectrum 1-Halo term. We looked 
at triangles for which k2 = 2fei — 20h~^ Mpc. Clearly, the 
analytic predictions and direct measurements from the syn- 
thetic data are in excellent agreement (2a level). We are 
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therefore confident that, for these problems, the integrators 
are accurate to at least ~ 1%. 

To summarize the main results of Section |^ we have 
convincingly demonstrated that the detailed shape of Q{0) 
on small scales depends very sensitively on the shapes of 
dark matter haloes, and also their axis ratio PDF. Further- 
more, we propose that the discrepancies observed between 
the theoretical predictions for Q{0) and direct measurements 
from A'^-body simulations may be largely explained through 
the break down in the spherical halo approximation. 



5 GALAXY CLUSTERING 

We now turn our attention to the second of the questions 
posed in the introduction, and explore to what extent halo 
shapes impact on the galaxy clustering. Such information is 
of great importance if one wishes to use the halo model as 
a means for constraining the parameters of the HOD from 
observations. 



5.1 Galaxies in the halo model 

As a number of authors have shown, the halo model for- 
malism can be extended to successfully model the cluster- 
ing of galaxies (Seljak 2000; Peacock & Smith 2000; SSHJ; 
Berlind & Weinberg 2002, Yang et al. 2003, Cooray & Sheth 
2003). The methodology follows ideas first laid down in ear- 
lier works by Ncvman & S ^Si^iitti ^ n LP ^^ ^^ '-^ others, but which 
are nicely summarized in lPeeblej l|li)80ll . 

We note at the outset that one can write down the hi- 
erarchy of correlation functions directly from the mass clus- 
tering statistics by making use of the following simple sub- 
stitutions. Firstly, we make the transformation 



U{r,M,Si,e) ^Ug{r,M,Ei,e) 



Pair) 



(48) 



where Ug describes the average spatial distribution of galax- 
ies within each halo normalized by the expected number of 
galaxies. Note that a sub- or superscript g labels galaxy 
quantities. Secondly, instead of weighting the integrals over 
the halo mass function by mass, we weight by either the 
expected number of galaxies, or galaxy pairs, triplets etc., 
depending on whether we are considering the 3-, 2- or 1-Halo 
terms respectively. Thus we have 



n(M)M 
n{M)M^ 
n{M)M^ 

where 



n{M) (N\M) 
n{M) {N{N - 
n{M) {N{N - 



1){N - 2)|M) 



{N\M) = ^NP{N\M) , 

JV=1 

oo 

{N{N - l)|Af) = ^ N{N - 1) P{N\M) 



(49) 
(50) 
(51) 



(52) 



(53) 



In the above equations P{N\M) is the conditional halo oc- 
cupation probability that gives the probability of finding TV 
galaxies in a halo of mass M. Lastly, in order to correctly 
normalize the clustering we perform the transformation 

where fig is the mean number density of galaxies. 
5.2 Galaxy correlation functions 

Following the procedure described in the previous section, 
we can write down the two- and three-point correlation func- 
tions for the galaxies. Thus, for the two-point function, we 
find to be: 

= C + ; (56) 

gg" = , / dMd:>idad£n(M){N(N -l)\M)p(a\M) 

ni Stt^ J 

2 

X Y[{Ug{M,r,-^,a,£)} ; (57) 

i=l 

/2 
T\{dM,dxidsi.d£^n{M,) {N\M,) 

X p(a,|M0(7g(A/i,ri -xi,ai,£:i)}C2(l,2) . (58) 
Similarly, we find to be: 

Cg = Cf + Cg'' + C|'' ; (59) 

(™ = / dMdxdad£n{M) {N{N - l)(iV - 2)|M) 

Tig- STT / 



X p 



(a|M)J];{(7g(M,r,-x,a,£:)} ; 



(60) 



_3,g / Y[ ^'^^^ ■niMi)p{ai\Mi) 

X U{^^-r„M„a„£,)}{N{N -l)\Ah} {N\M2) 
X (7(xi-r3,Mi,ai,£:i)C2(l,2)-Kcyc ; (61) 



1 r ^ 

^-2)3 J YildM^ dx, dai d£i {N\M,) n{M,] 



(ng8 



X p(a,|MOC/(x, -r„M„a„£:0} ^(1,2,3) . (62) 



5.3 Power spectrum &: bispectrum 

We next write down the power spectrum and the bispec- 
trum. As for the analysis of the dark matter bispectrum, we 
define the useful functions W^, Wf^ and Wf23 by applying 
the transformations described in Section 15.11 to equations 
I39I I. 14311 and 144|l . Hence, the galaxy power spectrum can 
be written as: 



p plH I p2H _ 



(63) 



{N{N ~ 1){N - 2)|M) = ^ iV(iV-l)(iV-2) P{N\M) .(54) 



Pg" = I dM {N{N ~ 1)|M) n(M) Wf^iM) ; (64) 
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Pi ' r 

= -T n h^^' {N\M,) n{M,) b{M,) W^iM,) . (65) 



Similarly, the galaxy bispectrum can be written as: 

ng _ nglH ng2H og2H . 
^123 ^ ^123 ^ -°123 ^ ^123 ' 

Bflf = ^ / dMn{M) {N{N ~ 1){N - 2)\M) 



(66) 



(67) 



Bf2f = ^Hj^ I n ^'^^^^ ~ l)!^i>} 



= {1,2} 



X {N\M2) lVf3(Mi)T^|(M2) +cyc ; 
^i^s"" = ;;3 / n^'^^' n(MOW^f(MO} 

/iff / 



B"i23{Mi,M2,M3) 



(68) 



(69) 



Owing to the fact that we are now modelling a discrete-point 
distribution, as opposed to a continuous field, we must apply 
corrections to both the power spectrum and bispectrum for 
shot noise. For this we follow SSHJ and use relations that 
have the same form as the Poisson case jPeeblejllQSol) : 



(70) 

cl23 — -"123 ~ '/L-" cl J c2 T c3J ~ '/ ' C^^) 

but where = P^(k oo) instead of the usual rj — l/fig 



^„1TJ — ^123 '/[-■ cl ' ^c2 ^ ^c3\ 'I 



pairs with double counting. The galaxy window function is 
then obtained by averaging the weight factor W over the 
probability distribution for the axis ratios and the Euler 
angles: e.g. in general the n-point window function is given 
by 



W?2..AM) = ^ j dadf p(a|Af)Wf2...„(M,a,f) . (73) 

On computing the averages over P[N\M) in equation 17211 . 
we obtain 

Wi^|(Af,a,f ) = {N{N - 1)\M) \U,{\^t 

+ 2 [(TV - 1|A/) + P(0|A/)] {t/,(k) - |t/3(k)|^} , (74) 

where we have suppressed the dependence of ?7g on a and 
£. The term P(0|Af), which is the probability of finding 
no galaxies in the halo, arises because the sum over galaxy 
numbers in equation 1721 does not start from zero. Hence 
we have terms like 



^(Af- 1)P{N\M) = {N - 1\M) -P{Q\M) 



(75) 



On combining equations 163ll - 165j . 1731 and 17411 . we arrive 
at a complete description of the galaxy power spectrum. 

Similarly, for the bispectrum we will have (A*' — 1)(A'' — 
2)/2 triplets that include the central galaxy and (A'^ — 1)(A^ — 
2)(A'^ — 3)/3! triplets that comprise only satellite galaxies. 
The 2-Halo term in the bispectrum will be modified exactly 
as described for the power spectrum. However, for the 1-Halo 
term the weight factor must be modified as follows: 



5.4 Central galaxy contribution 

The above equations assume that the spatial distribution 
of the galaxies follows the underlying dark matter profile. 
We may also wish to impose the further constraint that, for 
every halo that contains at least one galaxy, a single galaxy 
lies at the halo centre of mass. In terms of the halo model 
this means that there will be a small modification to the 
clustering for terms that include the correlation of objects 
within the same halo. We now derive how this affects the 
small scale clustering. 

Consider a triaxial halo of given mass A/, with axis ra- 
tios a, orientation £ and with A'^ galaxies. For this halo we 
can identify one central object and A'' — 1 satellites. For the 
two-point statistics, there will be (A'^ — 1)(A'^ — 2)/2 pairs of 
satellite galaxies and A" — 1 pairs that include the central 
galaxy and satellites. The two-point clustering of the satel- 
lite galaxies can be calculated as per usual. However, for the 
pairs that include the central galaxy, the clustering will be 
weighted by the density profile. This owes to the fact that 
the probability of finding a galaxy at position vector r from 
the central object now simply follows the density profile of 
galaxies. Hence, for the 1-Halo term in the power spectrum, 
we have the weight factor 

>V^|(Af,a,f ) = 2! ^ P(iV|A^) {(iV - l)(/,(ki, Af, a, £) 
Ar>l 

+ i(iV-l)(Af-2)t/,(ki,A./,a,f)C/;(k2,Af,a,f)| .(72) 

Note that the factor of 2! in the above expression arises 
since the term UgUg is normalized by the total numbers of 



Wj*|3(A/, a, f ) = 3! ^ P{N\M) i i(7V - l)(Af - 2) J] [/g(k.O 



-I 

+ -(iV-l)(iV-2)(iV-3)J][/g(kO 



(76) 



where we have included the factor of 3! to account for the 
fact that the profiles are normalized by the total triplets, 
which are counted in a way that does not prevent repetition. 
On computing the average over P{N\M) we get 



3(A//,a,£) = {N{N - 1){N ~ 2)iA/) [] t/g(k. 



+ 



3! 1 1 (iV(iV ~ 1)) -{N- 1\M) - P(0|Af)} 



(77) 



Thus on combining equations 1661 - 1691 . 1731 and 1771 we 

also arrive at a complete analytic description of the galaxy 
bispectrum with one central galaxy and A' — 1 satellites. 

5.5 Halo Occupation Distribution 

To calculate the galaxy bispectrum we require a specific 
model for the mean spatial distribution of galaxies within 
each halo and also the first three moments of the halo occu- 
pation probability function. Firstly, we will make the usual 
assumption that the galaxy density distribution follows the 
mass profile, and hence Ug{r, M) oc U{r, M). In accordance 
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Figure 8. Comparison of tlie first tliree moments of the halo 
occupation probability function of K04, solid lines, with those of 
SSHJ, dashed lines. In the top panel, the dot-dash line represents 
the central galaxy contribution and the dotted line represents the 
satellite galaxy contribution. 



with the ideas of the previous section, we wiU consider the ef- 
fects on the bispectrum caused by placing one galaxy at the 
centre of each halo. Secondly, for P{N\M), we will consider 
two potentially equally viable models from the literature: 
the first is the binomial model proposed by SSHJ, that was 
constrain ed to match the predict ions of the semi-analytic 
models of iKauffmann et al.l 1 1999|k the second is the model 
developed by Kravtsov et al 12004 . hereafter K04). In this 
model it is assumed that if the halo mass lies above some 
minimum mass threshold, then there is always one central 
galaxy. The remaining A'' — 1 galaxies are then taken to be 
satellites of the central galaxy, and it was found that these 



follow a Poisson process ^. The full details of both models 
are summarized in Appendix [HI 

In Figure |H| we show the differences between the first 
three factorial moments of the two probability functions. 
Note that in the limit of high masses, the power-law index 
of the first moment for the model of K04 is (N) oc M" , 
with a ^ 1.0, whereas for the binomial model it is a ~ 0.9. 
In this limit the second and third factorial moments scale 
identically with the mean (A'^) for both models. For halo 
masses below ~ 10^^ Mq (i.e. group and galaxy masses) the 
second and third moments are sub-Poisson in both models. 
Furthermore, for a given mass the binomial P{N\M) is much 
narrower than for the K04 model. 



5.6 Results: Galcixy Bispectrum 

Figures lUl and ITUl show the triaxial and spherical halo model 
predictions for Q{0) for the binomial model of SSHJ and 
the Poisson satellite model of K04, respectively. In generat- 
ing these predictions we chose the minimum halo mass for 
the binomial model to be Mmin = 10^^ Mq, which for 
the assumed cosmology gave a number density of galaxies 
fig = 3.81475 X 10~^fe^Mpc~^. In order to compare the two 
occupation probability models consistently, it was important 
to remove any scaling due to the overall numbers of galaxies, 
since Q oc fig . This was done by finding the Mmin in the K04 
model that reproduced the same galaxy number density as 
for SSHJ. We found this to be M^in = 1.28 x 10^^ h'^ Mq. 

The top panels of both figures show the predictions for 
the case where all galaxies are distributed within the haloes 
according to the density profile (Section 15.311 . The bottom 
panels show how the predictions change when one galaxy is 
assumed to be located at the centre of each halo (Section 
15.41 . As in Section |1| we show results for both the triaxial 
NFW profile (left panels) and the JS02 model (right panels). 
Once again, it is seen that the triaxial halo model predictions 
for Q{9) differ quite substantially from those of the spherical 
model. The characteristic features of these predictions are 
essentially the same for all of the eight panels in the figures: 
the triaxial model produces a boosted signal for colinear 
triangle configurations and a suppression for configurations 
that are close to equilateral. However, and importantly, the 
functional form of these curves appears not to be strongly 
dependent on the precise form of the chosen HOD. 

On comparison of these results with the mass bispec- 
trum (Fig. 2J , we see that the overall effect of triaxiality is 
weaker in the galaxy statistics. For the binomial model the 
effects are reduced by ~ 50%, and for the model of K04 they 
are reduced by roughly ~ 30%. The better sensitivity of the 
K04 model can be attributed to the fact that, for high mass 
haloes, the power-law index of the first moment of P(A'^|Af) 
is a steeper function of mass than the binomial model: as 
mentioned previously, a 1 for the K04 model and ~ 0.9 
for the binomial model. In other words the most massive 



^ Note that the term "central galaxy" used in K04 is somewhat 
unfortunate given our discussion in the last section, since the 
occupation probability model of K04 in fact makes no prediction 
for the spatial distribution of galaxies within the halo. Thus this 
"central galaxy" may, in principle, lie at any point within the halo 
boundary. 
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Figure 9. Configuration dependence of the galaxy bispectrum. The halo occupation probability function is taken to be the binomial 
model of SSHJ. Line styles are the same as in Figure HI 



haloes, which are also the most triaxial, will host relatively 
more galaxies in the K04 model than in the binomial model. 

Comparison of the upper and lower panels of Figures]^ 
and I lUI reveals the overall importance of having one central 
galaxy. The first point to note is that, for both P{N\M) 
models and density profiles considered, the overall effect is 
that Q{9) is shifted by a factor ~ 1. After some tests, we 
found that this behaviour was a consequence of the bispec- 
trum being normalized by products of the power spectrum 
("equation I47I I. To see this consider the following: on small 
scales, the power spectrum is more strongly weighted to con- 
tributions from low mass, low occupancy halos; this is due 
to the {N{N - 1)|M) weighting in the 1-Halo term. Whilst 
both the power spectrum and the bispectrum are sharply 
enhanced by the presence of the central galaxy, the power 
spectrum is relatively more enhanced owing to the added 
numbers of lower mass haloes that contribute to the signal; 
these are not included in the bispectrum as only haloes that 
host at least three galaxies can contribute. The result then 
follows. 



Secondly, the central galaxy also appears to reduce the 
overall affect of halo triaxiality on the clustering. This can be 
understood by considering haloes that host three galaxies: 
when one is placed at the centre we find that the 1-Halo term 
is B(k) oc [Ug{k)]^ (equation |7SJ, whereas for the case of no 
central galaxy, B(k) oc [Ug{k)]^. As was shown in Figure 
the window functions become more sensitive to the halo 
shape as more points are placed in the fc-space halo. Thus, 
placing one galaxy at the halo centre effectively removes one 
of these points. 

As a final remark, we note that the combination of mea- 
surements for the triangle configurations with ^2/^1 = 1, 2 
and 3, in both the spherical and triaxial models, appear to 
depend strongly on the HOD. In particular, the spread in the 
amplitude of Q as the ratio ^2/^1 changes is substantially 
greater for the model of KS04, while the predictions for the 
SSHJ model appear closer to those of the mass bispectrum 
(cf Figure^. The response of these predictions to the pres- 
ence, or otherwise, of a central galaxy is also very different 
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Figure 10. Same as Figure|2l however, this time the halo occupation probabiUty function is taken to be the Poisson satellite model of 
K04. Again, line styles are the same as in Figure El 



in the two models, with the SSHJ model being aftected to a 
much greater extent by this assumption. 

In summary, on small scales the reduced galaxy bis- 
pectrum has been shown to be a sensitive function of the 
halo shape and density structure, the occupation probabil- 
ity P{N\M), and the central galaxy assumption. We note 
that the functional form of Q{6) is most strongly affected 
by the shapes of the haloes themselves and is not greatly 
affected by the exact form of the HOD. However, the am- 
plitude of each prediction, relative amplitudes of predictions 
made over different scales and the strength of the triaxiality 
effect are all sensitive to the particular choice of HOD. 



6 CONCLUSIONS 

In this paper we have explored how sensitive the real space 
bispectrum in cosmology is to the underlying shapes of dark 
matter haloes, for both matter fluctuations and galaxies. We 
achieved this through application of the triaxial halo model 
of SW05, adapting it to the problem of higher order statistics 



in a straightforward manner. In our model we accounted for 
dark matter haloes that have some known distribution of 
shapes and that are randomly orientated in space. 

Analytic expressions were written down for the real 
space three-point matter correlation function and bispec- 
trum. We then concentrated on the bispectrum and found 
that the 1-, 2- and 3-Halo terms could be compactly ex- 
pressed as 7-D integral equations. These were solved numeri- 
cally using an efhcient multi-dimonsional quadrature routine 
(Korobov 1963; Conroy 1967). Two models for the density 
profile of the triaxial haloes were considered: the first was 
the triaxial NFW model of SW05, which allowed us to ex- 
plore the affects of halo shape alone on the bispectrum; the 
second was the more realistic model of JS02. The bispectrum 
was examined over a wide range of scales and for different 
configurations of fc-space triangles. 

We found that for equilateral triangle configurations the 
effect of halo triaxiality on the bispectrum was to produce 
a suppression on scales k > 0.2/iMpc~^, relative to the 
spherical halo model predictions. This suppression was at 
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the level of < 7% for both density profiles. The reduced 
bispectrum was considered next, as function of fc-space tri- 
angle configuration. We found that on large scales this quan- 
tity was insensitive to halo shapes. However, on small scales 
the predictions for the spherical and triaxial halo models 
differed significantly in both the amplitude and functional 
form of Q{9). In general we found that, relative to the spher- 
ical model, the triaxial halo model predictions produced ex- 
cess signal for triangle configurations that were colinear and 
non-vanishing, and a deficit for triangles that were close to 
equilateral. We observed that the overall form of Q{9) de- 
pends little on the density profile model itself, modulo an 
amplitude shift, but rather that it is the shape of the haloes 
that dictates the shape of the curves. 

Our results were then considered in relation to the dis- 
crepancies between the A'"-body simulations and the predic- 
tion s from the spherical halo model that were noted by SSHJ 
andlPo salba. Pan fc Szaoudil (I2flfl 5^ . We suggested that, for 
the case of the equilateral triangle configurations, the in- 
consistencies between theory and simulations are not due 
entirely to a break down in the spherical approximation. 
Rather, they are likely due to the combination of finite vol- 
ume effects in the simulations and also other corrections to 
the halo model (SSHJ; Wang et al. 2004; Fosalba, Pan & 
Szapudi 2005). However, for the case of the configuration 
dependent bispectrum it was shown that the discrepancies 
could be well attributed to the break down of the spheri- 
cal halo approximation, modulo an amplitude offset, which 
arises from the correction required to reconcile halo model 
predictions for the equilateral bispectrum with numerical 
simulations. 

In the second part of this paper, we explored how halo 
shapes affect the real space galaxy bispectrum. Galaxies 
were included in the halo model in the usual way (Seljak 
2000; SSHJ; Berlind & Weinberg 2002), and we considered 
two different schemes for the halo occupation probability 
function. The first of these was the binomial model devel- 
oped by SSHJ, while the second was the Poisson satellite 
model of K04. For both probability functions, the predic- 
tions of the configuration dependent, reduced bispectrum 
were found to be very similar to those for the mass. However, 
the overall magnitude of the effect was reduced by ~ 50% 
for the binomial model, and ~ 30% for the K04 model. It 
was argued that the K04 model was more sensitive to halo 
shapes because of the stronger dependence on halo mass 
for the first moment of the occupation probability function. 
We also showed that the predictions for the binomial model 
possessed a stronger dependence on the presence of a central 
galaxy, when compared to the model of K04. Importantly, we 
found that the essential effects of triaxiality on the reduced 
bispectrum did not depend strongly on the exact details of 
the halo occupation distribution. However, a combination 
of measurements on different scales are sensitive to the both 
halo shape and the HOD. 

We conclude that the bispectrum provides significant 
information about the shapes of dark matter haloes in ad- 
dition to the (galaxy and mass) halo density profiles, and 
the halo occupation probability function. Consequently, in 
order to use the halo model to make precise constraints on 
the HOD, it will be necessary to take into careful considera- 
tion halo triaxiality. We remark, however, that our study did 
not account for the effect of redshift space distortions, which 



will have a significant impact on the behaviour of the galaxy 
bispectrum, particularly on small scales. Furthermore, we 
have also neglected the changes to the bispectrum that are 
caused by selecting different galaxy populations i.e., popu- 
lations based on colour, luminosity or type. Neither have we 
considered realistic selection functions, or survey geometry. 
We reserve such investigations for future work, when we will 
study in more detail the bispectrum as a tool for measuring 
the HOD. 

It remains to be seen whether measurements of the 
three-point statistics in weak lensing will be sensitive to halo 
sh apes. An initial stud y of this problem has been conducted 
bv lHo fc Wliit3 ll2004l) ■ who calculated the three-point shear 
correlation function for triaxial haloes and found strong ef- 
fects. The situation they considered was rather idealized, 
however, and was not a full halo-model calculation. Nev- 
ertheless, the p otential for mea s uring s uch effects is clear. 
Recent work bv iTakada fc Jaiiil ll2003tfl showed, for an 11 
square degree mock cosmic shear survey, that the three- 
point shear correlation function from ray-tracing simulations 
was in good agreement with the predictions of the spheri- 
cal halo model. However, they noted that on small scales 
there were significant and unexplained discrepancies. For 
forthcoming weak lensing surveys which aim to cover even 
larger fractions of sky (> 100 square degrees), the shear 
three-point functions will be accurately determined over a 
wide range of scales. Discrepancies between halo model pre- 
dictions and observation will become very significant and 
important to understand. Moreover, the detection of shape 
information in the mass bispectrum or correlation function 
via weak lensing would be an important consistency test of 
the CDM paradigm, since it has been shown that plausible 
variants of the dark matter model, like SIDM and WDM, 
will giv e rise to a different spectrum of dark matt er halo 
shapes jAvila-Reese et alJl20oil lYoshida et al.ll2000ll . 

Lastly, we mention that a recent study of the shapes of 
da rk matter haloes in numer ical simulations with gas cooling 
bv iKazantzidis et alj (|200^ has shown that the formation 
of central condensations may act to sphericalize the inner 
regions of haloes. At the virial radius, where halo shapes are 
usually measured, such processes will have relatively little 
effect. However, in order to achieve precise predictions to 
compare with observation, it will be necessary to quantify 
in detail how these effects modify halo profiles and axis ratio 
distributions. 
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APPENDIX A: ROTATION MATRIX 

Owing to there being several equivalent ways to define the 
Euler angles for the rotation matrix TZ{a, 13,"/), we make 
explicit the definition that we use throughout. The matrix 
for the z ~ y' — z" rotation is given by: 



7^(a,A7) = 



{ Cp Co. 

Sol C'y^ 



(C*/3 Sa C^ 
H-Ca 5''y) 
( Cp Sa. S^ 
-\-Ca C'-y) 
Sq Sa 



— Sp C-y 
Sp S-y 



.(Al) 



where we have adopted the short hand notation Cx 
and Sx ~ sin s. 



APPENDIX B: OCCUPATION PROBABILITY 

Bl Scoccimarro et al: binomial distribution 

We now summarize the details of the occupation probabil- 
ity model proposed by SSHJ to describe t he results from the 
semi-analytic galaxy formation models of iKauffmann et alJ 
l|l999^ . ISheth fc Diafericl l|200ll ) measured the first two mo- 
ments of P{N\M) from these models. For the first moment 
of they found 



(iV|A/> = {Nb\M} + {Nr\M} 
M 



(Bl) 
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where Nb and Na represent the numbers of red and blue 
galaxies. The blue galaxy parameters are q_b = for haloes 
with masses in the range IO^^Mq/i"^ < M < Mb and as = 
0.8 for M > Mb, where Mb = 4.0 x IO^^Mq/i'V The red 
galaxy parameters are an = 0.9 and Mb = 2.5 x IO^^Mq/i"^ 
for M > W^^ Mq. For the second moment, they found 
that a sub-Poissonian distribution was preferred for low- 
mass haloes and an almost Poisson distribution for higher- 
mass haloes. They characterized this through 



{N{N ~ 1)\M) = A^(M) {N\Mf 



(B3) 



where the function A expresses the deviation from a Poisson 
probability distribution and has the form 



A(M) 



logio(M/Mii)'/' 
1 



M < 10" Mq h- 
M > 10"Mq h- 



(B4) 



where Mn = lO^M© h'K 

SSHJ then proposed that, since P{N\M) was not well 
described by a Poisson process, a better description might 
be afforded through the binomial distribution: 



P(N^n\M)^C^'''p"Mil~PM) 



(B5) 



Nm\/{Nm — n)\n\. The binomial probabil- 



where Cn"^^ 

ity function is completely specified by two free parame- 
ters: the probability that out of one trial the event occurs 
Pm = p{M); and the total number of trials Nm = iV™'"'(M). 
These parameters were then fixed by matching the first and 
second moments to those measured from the semi-analytic 
galaxies. This gives 



Pm = {N\M) \l - \\M) 



Having specified P(A''|A/) one may then calculate the 
galaxy clustering to any desired order. Of particular inter- 
est are the factorial moments. For the binomial distribution 
these are most readily obtained by differentiation of its fre- 
quency generating function: 



{N{N ^1)...{N ~j)\M) 



d'Fit) 



dp 



= j^iPMi + l-pM}""" 



(B7) 



On performing the differentiation and substituting in for the 
relations <B6ll . one finds the general relation (SSHJ) 

{N{N ~1)...{N -j)\M) = 

A"(2A2 - 1) . . . (jA' - j + 1) {N\My+^ . (B8) 

As a final remark we note that this particular occupation 
function is somewhat problematic to implement in practice: 
in order to sample from it, one requires that Nm be an 
integer. In an analytic calculation this is not the case since 
one obtains the factorial moments directly from equation 
ijBSfl . Indeed, as A is a continuous variable, so to will be 
Nm ■ Thus when Nm is small there are significant differences 
between the sampled and analytic distributions. 



properties of these objects as distinct populations (K04). 
The occupation probability under this prescription is there- 
fore re-written 



P{N\M) = J2 Pccnt{l\M)P,^t{N - l\M) 



(B9) 



where Pcant{l\M) and Psa.t{N — l\M) are the probabilities 
for getting I central galaxies and N — I satellite galaxies, 
respectively. 

In the work of K04 it was supposed that there was a 
one-to-one mapping between haloes (parent or substructure) 
and galaxies (central or satellite), and also that there was 
always only one central galaxy associated with the parent 
halo. This then lead to the following statement: 



Pccnt(iV = 1\M) = 







M > Mmin 
M < Mmin. 



(BIO) 



Under this condition equation l|B9|l simply becomes 

P{N + 1\M) = Peat(iV) ; M > Mmin ■ (Bll) 

Following this one may write down how the moments of the 
full halo occupation probability P{N\M) are related to the 
moments of the satellite population Psat(A''|M). 

K04 went on to show that the first moment of the satel- 
lite galaxy probability function was well characterized as a 
single power-law, {N)^^^ cx M, with a steeper decline for 
lower halo masses. The exact function they found was 

(iVU. = (^-c)\ (B12) 

where /? = 1.03, C « 0.045 or Mi/Mmin = 22, where M^in 
is the minimum halo mass. Furthermore, it was found that 
the second and third moments of the satellite galaxy distri- 
bution were well characterized by a Poisson process. Thus, 
the first three factorial moments of the full halo occupation 
probability function are: 



(N) = (iV,) + 1 ; 
{N{N ~ 1)) = (TV,)" 
{N(N -\)(N -2)) 



{Ns 



3 (Ns) + 2 



(B13) 
(B14) 
(B15) 



B2 Kravtsov et al: central and satellite split 

Much recent attention has focused on the separation of cen- 
tral and satellite galaxies, and characterizing the statistical 
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